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1. INTRODUCTION 


Consider the optimization problem ^ 

min f(Xj , . . . , Xj^ ) ; subject to ^ x. = 0, 

(1) i=l 

X. S 0 , i = l , . . . , n 

Problems of the form (1) occur frequently in statistical applications. 
In " se applications, the variables x^ usually represent a discrete 
probability density function and the functional f represents an 
optimality criterion, e. g. , the negative likelihood or negative log 
likelihood. 

Let us make the change of variables 

(2) X. = y? , i=l , . . . , n 


and consider the optimization problem 


n 


where 


(3) min f (y , . . . , y ) *, subject to 1 - ^ % yf = 0 

i=l 


<4) t(yi y„) = f<iyi iy„) 


Proposition 1. 1 

T 

y = (y » . . . » yj^) solves problem (3), then 
2 2 T 

x = (iy^,..., iyjj) solves problem (1). Conversely, if 
X = (x , . . . , x^) solves problem (1), then y = (/25c^ > • . . , 
v/ 2x^)"^ solves problem (3). 

Proof 


The proof is not difficult. 


} 


1 


I 



2 




Problem (3), in contrast to problem (1), does not 
involve inequality constraints. The price one pays for this simpli- 
fication is that problem (3) will have nore critical points than problem 
(1) and the equality constraint is quadratic instead of linear. We 
maintain, that according to Proposition 1.1, if one starts sufficiently 
near the solution, the extraneous critical points will not affect the 
algorithm. Moreover, when f is nonlinear, the quadratic constraints 
should not significantly increase the degree of complexity of the problem, 
Tapia’s approach to quasi- Newton methods for constrained 
optimization ([l] and [2]) depends heavily on the invertibility of the 
Hessian of the Lagrangian at the solution. Hence, we are certainly 
interested in determing the relationship between the invertibility of 
the Hessian of the Lagrangian for problem (1) and the invertibility of 

the Hessian of the Lagrangian for problem (3). Toward this end let 

n 

' (5) L (x,u,X) = f(x) - Y, ^^i ^i "I ^i^ 

i= 1 


and 

A n 

(6) L (y,^) = f(y) + X(1 - ^ I yf ) 

i = l 


we use the notation 


and 


A 

Y 



diag (y 1 , . . . ,7^) 


diag (x^ , . , . , x^) 


X 


Straightforward calculations show that 


(7) Vy L(y, X) = Y (vf(x) - a) , 

(8) L(y, X) = Y v^f(x)Y + diag (v f(x) - A 

y 

(9) L(x, u, X) = V f(x) - A - u , 
and 

(10) L(x, u, X) = f(x) . 

*ne first order necessary conditions for problem (1) 
are 

vf(x)-'A-p = 0 

1 - 5^=0 
i = l 

( 11 ) | i .&0 

X . 0 

l-tjXj^ — 0 I i — ly.ti^n f 

while the first order necessary conditions for problem (3) are 

Y ( V f(x) - A ) = 0 

(12) n 

1 - E iy? =0 

i = l 

Recall that strict complementarity means that in (11) we do not 
have li. = Xj = 0 for any i , 
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Proposition 1. 2 

Let (x, u, X) be a solution of problem (1) and (y, X) 
the corresponding solution of problem (3), Then 

2 

(i) we have strict complementarity and v f(x) is 

o A 

positive definite <==>v" L(y, X) is positive definite 

^ 2 

(ii) we have strict complementarity and v f(x) is 

o 

invertible <=> in' y L(y, X) is invertible. 

Proof 

Again the proof is not difficult. One merely works with 
(8), (10), (11) and (12). 

Proposition L 2 is very satisfying and warns us tha in 
the algorithm we must guard against the loss of strict complemen- 
tarity. This will be accomplished by working with the variable 


(13) 



where E. = 2 if jx. [ + (U J 9 ^ 0 and E, = 1 if 
[Xj^l + [Ujl =0 . Observe that, when E. = 1 in (13), we 
have effectively ignored the nonnegativity constraint on x^ . 

II. THE QUASI-NEWTON ALGORITHM 

The algorithm we are about to present is a special case of 
the approach given by Tapia in [2] applied to problem (3). The 
reader is referred to that paper for a better understanding of the 
theory and algorithm. In describing an iterative method, it is 
convenient to denote the present iterate by x and the subsequent 
iterate by x and similarly for other variables. Also, we let I 


f 


5 
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denote the identity matrix , •) denotes the Euclidean inner 

product, =: (1,0, . . . > ' ^ e^ = (0 , . . . , 1 )^, 

(a^j) denotes the matrix whose (i,j)-th element is a^^j 
and or (b). denotes the i-th element of the vector b . 
Consider the constrained optimization problem 

(14) min f(x), subject to g(x) = 0 ; 

where f , g : _* , We first describe the quasi- Newton 

method BFGS given in [2] for problem (14). Let 


(15) L(x, \) f(x) + \ g(x) 


Quasi- Newton BFGS 

Step 1 : (Initialization) Determine x and B . 

Step 2 ; (Update x and \) Let 

— = g(^) - < w g(x), 7f(x) ) 

( tr g(x), B" ^ V g(x) > 

X = X - B“^ L(x, 1 ) . 


Step 3 ; (Update approximate Hessian) 


where 


and 


Step 4; 


B = B + 


y y 


T 


< y » s > 


B s s'^ B 
(s, Bs ) 


\ 

I 


s = X - X 


y = L(x, x) - 7 ^ L(x, \) . 
GO TO STEP 2 . 



In implementing this algorithm, one would have to make pro- 
visions for output and a stopping criterion. In [2], Tapia has demon- 
strated that the above algorithm is locally superlinearly convergent. 
Remark 1 ; 

If the initial B depended on \ , we would 
m'Ske an initial approximation of x according to the projection 
formula (see [2]’ 

(16^ X = - f ^ g(^) ■ ^ > 

(V g(x) . 7 gOc) > 

Remark 2 : 

In updating B in Step 3, it is advantageous to 
store B factored according to a Cholesky decomposition and then 
merely update these factors. Moreover, B can be slightly 
modified so that the resultant matrix is always positive definite. 

For a detailed description of these modifications, the reader is 
referred to Gill et al [3], Van Rooy et al [4], Fletcher and 
Powell [5] and Van Rooy [6], The quantities B'^ ^ f(x) and 
B"^ ^g(x) in Step 2 are then obtained in the standard back- 
substitution manner. 

We now apply the above algorithm for problem (3) with the 
additional feature of obtaining a discrete Newton approximation to 
the Hessian every M iterations (where M is an input parameter). 

Define the following functions 


( 17 ) 


1 


I 

I 


(18) 


= - [2 - Ej + (E; - 1) yj ] , 




(19) g = 1 - X 

i = l 

The quantity y g is merely the gradient of g with respect to y 
and (13) and (17) are the same. 

Quasi" Newton Algorithm for Problem (3) 

Step 1 ; (Initialization) 

Read CRC (constraint regection criterion) 

h (discretization step) 

e (stopping criterion) 

DI (descent indicator) 

I0UTPUT (output indicator) 

M (discrete Hessian indicator) 

MAX (maximum number of iterations) 


E. = 2 


X = ( 7 f(x) , X > 

B = I 


Step 2 : Output according to iQUTPUT, 

check stopping criterion 



1 


I 

t 



» 






o 


¥ 


8 


IP (ll-’y L(y. >>) II + < e or 

ITEKATION > MAX ) stop 


Step 3 ; 


Calculate discrete Hessian according to M 


B = 



A 

L ( y + ® j * 



- V Uy, 



Calculate modified Cholesky decomposition of B . 
Step 4 : Calculate multiplier correction 

g - < <7g , B" ^ L(y , x) > 


< V g , B"^ yg > 


S tep 5 : 


Step 6 : 


Calculate y-variable correction 

Ly- - Vy L(y, x) - A x B"^ n g 

Update y according to descent indicator 
X = X + O’ A X 
y = y + O' A y 

where a - I if descent is not desired, other- 
wise 01 is chosen according to the descent 
principle on || Vy L (y , - x ) || ^ + 

(see [2]). 


Step 7 : Update approximate Hessian according to BFGS 

and modified Cholesky decomposition 



Nl 


i 

1 


Step 8 : 


h 


Step 9; 


B = B + Py - BPS Ds B 

(Dy,Ds > (Ds>BDs> 


where D s = x - x and 


P y = Vy L(y, X) •' L(y , x) 


Update constrain*' rejection indicators 
Tq = [ll Vy L(y, X) !! + 

(1- 1 ='-■ r ] 

i = l 


2.^2 


Ti = (?¥■ '^0 ) - 


T 2 = min (CKC , ) 


If (J Xj I < T and 


3 f 


I 3 X. " ^ I < ^^2 ) ' 
1 ' 


let E . = 1 


If has changed, reset y. 


i = 1 , . . . , n 


GO TO STEP 2 . 


a 


I 


{ 
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III. USER'S GUIDE 


The above algorithm has been implemented in a set of 
FORTRAN subroutines whose entry point is the subroutine CONOPT, 

In this section, we shall describe how to use this program, provide 
a description of the workings of the program, ard illustrate its use 
with examples. A listing of the program is given in the appendix. 

This program runs on an IBM 370/155 using the FORTRAN G , 

G1 , or H compilers. All computations are done in double precision. 
Unit 6 is used for output, if requested. 

To utilize this program, the user supplies several parameters 
and two subprograms which calculate the function value and its gradient 
at a specified point. The program will then attempt to minimize the 
function subject to the constraints and return the argument at the 
minimum. 

The calling sequence of the program is CALL CONOPT 
(N, DELF, FUNCT, X, WKA , MITER, STC , CRC, IDN, DH, 

KZJUT, DSNT, A, RSTRT, USELAM, XLAM) . 

These parameters are described below : 

N - integer variable indicating the dimension of the 

problem. 

DELF - a user supplied subroutine which calculates the 
gradient of the function. The calling sequence 
is CALL DELF (DF, X, N) , where X 
is the N dimensional argument and DELF 
returns the gradient in the vector DF . Note 
that DF and X are both double precision 
variables. 



FUNCT “ 


X - 




I 

I 
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a user supplied double precision function sub- 
program which calculates the function value at 
a specified point. Tlie cplling sequence is 
F = FUNCT (X, N) where X Is the (double 
precision) argument of length N . Though the 
function value is not required by the algorithm, 
CONOPT will output it in the convergence test 
step if the user requests output. (It is called 
only when output is requested.) If, for any 
reason, the user is not interested in this value, 
he may just supply a function subprogram that 
sets the function value to zero and returns, 

[N. B. The two programs DELF and FUNCT 
must be declared external in the calling program. 
The two programs may have other entry points, 
with different arguments, referred to by the 
user's calling program prior to his calling 
CONOPT in order to make other variables or 
their addresses available to either or these 
subprograms. ] 

on output, this will contain the N- dimensional 
double precision argument last calculated by the 
program (at the minimum if the program 
converged). On input, the program will set 
each element of X to the value 1/N , unless 
the user has specified restart (see description 
of the variable RSTRT below), in which case 


12 





I I i 


* 


K 


WKA - 


MITER - 


STG - 


CRC - 


IDN - 


the supplied values of X will be used as the 
starting point. 

a double precision working storage array of 

length &^N^ + (23N + 7)/2. If a restart 

has been specified, the first N (2-byte integer) 

locations should contain the desired constraint 

rejection indicators. However, if these variables 

do not contain either I's or 2's , the program 

will set them to 2. The final constraint rejection 

indictitors will be in these locations. 

an integer variable indicating the maximum number 

of iterations to perform. If the user sets this to 

0, the program will reset it to max (30, 2 * N). 

a real variable specifying the stopping criterion on 

the norm squared of the gradient of the Lagrangian 

plus the constraint, i. e. || ^ L ((^ + (l " T ^i) "• 

i = l _o 

If 0 is specified, the program will set STG = 10 . 

a real variable specifyi the constraint rejection 
criterion. See the algorithm description section 

for a precise definition of its use. If 0, is 

_2 

specified, the program will set GRG = 10 
an integer variable specifying the frequency of 
calculation of the discrete approximation to the 
Hessian, i. e., this finite difference Hessian is 
calculated every IDN iterations. If 0 is 
specified, this approximation is never used. 
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DH 


IC^UT - 


DSNT - 


A - 


If a negative value is specified, the program 
sets IDN = min ( 10 , N) . 
a real variable specifying the step size to be 
used in calculating the finite difference Hessian 
(see the algorithm description for a detailed 
definition). If 0 is specified, the program will 
set DH = 10 . 

an integer variable specifying the frequency of 
output ; L e. , the program will output Various 
computed quantities every IC)UT iterations. 

If 0 is specified, no output will be produced. 

If a negative value is specified, the program 
sets I0UT = 1 . 

a one- byte logical variable indicating whether the 
program should implement descent in its search 
at each iteration. 

a two-word real vector containing respectively 
the multiplier and exponent for DNORM 

(=|!vLl[^ + (l-^ ^i)^) 

i = l 

constraint rejection test ; i. e. , the appropriate 
quantity is A(l) * DNORM * * A(2) . The 
defaults for A(l) and A (2) are 1. and 
. 2 respectively. If the ijse.r wishes to set 
A ( 2 ) = 0 . , he should set A ( 1 ) to the 
negative of what he wants for A ( 1 ) and 
then a value of zero will be used for A ( 2 ) . 



RSTRT - a one- byte logical variable indicating whether the 
user is supplying an initial value for X (as in 
a restart case). If RSTRT ^ . TRUE , , the 
first 2N bytes of WKA should contain the 
vector ICR (see description of WKS and the 
algorithm description— the E/s ). 

USELAM - a one-byte logical variable used if RESTR = 

. TRUE, to indicate if the user is supplying a 
value for XLAM . If RESTR = . TRUE . and 
USELAM = . FALSE . , the value < v f , x > 
is used for XLAM . 

XLAM - the Langrange multiplier for the equality constraint. 

Used on input only if USELAM = RSTRT = . TRUE . 

The default values for the above parameters have boen 
chosen so that many problems will converge; thus, the only input 
parameters with which the user need concern himself in many cases 
is N, ie)UT, DSNT and RSTRT (usually =. FALSE . ) 

The subprograms used and a brief description of their function 
is given below : 

CONOPT - main driver subroutine 

DIMS - computes base addresses for various working 
arrays used throughout the calculations, 

INPT - sets default values for parameters and outputs the 
parameters used. 

INIT initializes y . x , ICR , H and x . Also fetches 

values of yf and vL . 



GHKSTP - 
DHESS - 

DEL A MY ■ 
NEWY - 

UPDTH - 

CONSTR ' 

GRADE - 
SUM - 

MCHLSK - 
PDSOLV - 


checks the stopping criterion and outputs quanti- 
ties if requested. 

calculates a finite difference approximation to 
the Hessian and its modified Cholesky decompo- 
sition, if requested. Tlie Hessian is forced to 
be positive definite, 

• computes the update terms v \ and n y , 

updates y and \ based optionally on descent 

of II VL 1!^ + ( 1 - ^ X. ) . 

i = l 

updates the approximate Hessian using the 
Broyden- Fletcher- Goldfarb- Shannon (BFGS) update 
technique. (Really the MCD of it is updated. ) 
Here again the Hessian is forced to be positive 
definite. 

updates the constraint rejection indicators and, 
if necessary, y and V L . 
computes v L given v f , x , y , and ICR . 
a function subprogram that performs various 

utility functions ( £ ’ 

V \ i = l i=l 

1 ^i ^i )' 
i = l 


computes the modified Cholesky decomposition 
(MCD) of a matrix. If the matrix is not 
positive definite, the MCD is modified to 
represent a matrix that is (see e. g. ref 4). 
solves Ax = y where A contains the MCD 


of a matrix. 


COMPT - 


computes a rank one update of an MOD usin^ 
f e composite- 1 method (see refs, 5 and 6). 
LDMUL - multiply a vector by a. matrix whose MCD is 
supplied. 


Examples ; 

The following quadratic programming (i. e. , min || A x - b H 
3. t* 0 and = 1 ) examples were run using all the 

default options ; 


1 . 


A = 



b = 



X = 


(.11112108, .44443959, .44443959) 


congerged in 5 iterations. Exact answer == 
(1/9, 4/9, 4/9)’’^ 


The robustness of the algorithm is exemplified by the 
following two examples for which the algorithm, theoretically, 
should not converge since the Hessian is singular at the solution. 


2 . 


A 




X = (.50000019, 0., .50000004)^ 


converged in 14 iterations. Exact answer = 

(%» 0 , 

The constraint for X 2 is not active in this case. 
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congerged in 11 iterations. The constraints for Xj^ 
and X 2 are not active in this case. 
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Program Listings of CONOPT 
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